function LEs = LE(a, b2, c, x0, y0, z0, N, transient, delta)
    X = [x0; y0; z0];
    Q = eye(3);
    sum_LE = zeros(3,1);

    for i = 1:N
        % 直接推进，无极小值保护
        X_next = SIMM(X(1), X(2), X(3), a, b2, c);

        % 数值差分法雅可比
        J = Jaco(X(1), X(2), X(3), a, b2, c, delta);

        % QR正交
        Z = J * Q;
        [Q, R] = qr(Z);

        r_diag = diag(R);  % 不做任何修正

        if i > transient
            sum_LE = sum_LE + log(abs(r_diag));
        end

        X = X_next;
    end

    LEs = (sum_LE / (N - transient))'; % 输出行向量
end

